<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of mvk</title>
  <meta name="keywords" content="mvk">
  <meta name="description" content="DMRG fast matrix-by-vector product, Method 1 (less accurate)">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html v1.5 &copy; 2003-2005 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../index.html">Home</a> &gt;  <a href="#">tt2</a> &gt; <a href="index.html">@tt_matrix</a> &gt; mvk.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../index.html"><img alt="<" border="0" src="../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="index.html">Index for tt2/@tt_matrix&nbsp;<img alt=">" border="0" src="../../right.png"></a></td></tr></table>-->

<h1>mvk
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>DMRG fast matrix-by-vector product, Method 1 (less accurate)</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>function [y]=mvk(a,x,eps,nswp,z,rmax) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre class="comment">DMRG fast matrix-by-vector product, Method 1 (less accurate)
   [Y]=MVK(A,X,EPS,NSWP,Y,RMAX) Matrix-by-vector product of a TT-matrix A
   by a TT-tensor X with accuracy EPS. Also, one can specify the number of
   sweeps NSWP, initial approximation Z and the maximal TT-rank RMAX (if
   they become too large)


 TT-Toolbox 2.2, 2009-2012

This is TT Toolbox, written by Ivan Oseledets et al.
Institute of Numerical Mathematics, Moscow, Russia
webpage: http://spring.inm.ras.ru/osel

For all questions, bugs and suggestions please mail
ivan.oseledets@gmail.com
---------------------------</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="../../tt2/@qtt_tucker/diag.html" class="code" title="function [qt]=diag(qt)">diag</a>	Diagonal of a matrix or diagonal matrix from a vector in QTT-Tucker</li><li><a href="../../tt2/@qtt_tucker/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>	Frobenius norm of the QTT-Tucker</li><li><a href="conj.html" class="code" title="function [b]=conj(a)">conj</a>	Complex conjugate of a TT-matrix</li><li><a href="diag.html" class="code" title="function [tt]=diag(tm)">diag</a>	Extract the diagonal of the TT-matrix</li><li><a href="norm.html" class="code" title="function [nrm] = norm(t,varargin)">norm</a>	Matrix norm of the TT-matrix</li><li><a href="size.html" class="code" title="function [sz] = size(tt)">size</a>	Mode sizes of the TT-matrix</li><li><a href="../../tt2/@tt_tensor/conj.html" class="code" title="function [tt1]=conj(tt)">conj</a>	Compute a complex conjugate of TT-tensor</li><li><a href="../../tt2/@tt_tensor/diag.html" class="code" title="function [tm]=diag(tt)">diag</a>	Constructs diagonal TT-matrix from TT-tensor</li><li><a href="../../tt2/@tt_tensor/norm.html" class="code" title="function [nrm] = norm(tt)">norm</a>	Frobenius norm of the TT-tensor</li><li><a href="../../tt2/@tt_tensor/qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>	Left and right orthogonalization of the TT-format</li><li><a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>	Reshape of the TT-tensor</li><li><a href="../../tt2/@tt_tensor/size.html" class="code" title="function [sz] = size(tt,dim)">size</a>	Mode sizes of the TT-tensor</li><li><a href="../../tt2/core/my_chop2.html" class="code" title="function [r] = my_chop2(sv,eps)">my_chop2</a>	Truncation by absolution precision in Frobenius norm</li></ul>
This function is called by:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="mtimes.html" class="code" title="function [c] = mtimes(a,b,varargin)">mtimes</a>	C=A*B: Matrix-by-matrix, matrix-by-vector, matrix-by-number multiplication</li></ul>
<!-- crossreference -->



<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [y]=mvk(a,x,eps,nswp,z,rmax)</a>
0002 <span class="comment">%DMRG fast matrix-by-vector product, Method 1 (less accurate)</span>
0003 <span class="comment">%   [Y]=MVK(A,X,EPS,NSWP,Y,RMAX) Matrix-by-vector product of a TT-matrix A</span>
0004 <span class="comment">%   by a TT-tensor X with accuracy EPS. Also, one can specify the number of</span>
0005 <span class="comment">%   sweeps NSWP, initial approximation Z and the maximal TT-rank RMAX (if</span>
0006 <span class="comment">%   they become too large)</span>
0007 <span class="comment">%</span>
0008 <span class="comment">%</span>
0009 <span class="comment">% TT-Toolbox 2.2, 2009-2012</span>
0010 <span class="comment">%</span>
0011 <span class="comment">%This is TT Toolbox, written by Ivan Oseledets et al.</span>
0012 <span class="comment">%Institute of Numerical Mathematics, Moscow, Russia</span>
0013 <span class="comment">%webpage: http://spring.inm.ras.ru/osel</span>
0014 <span class="comment">%</span>
0015 <span class="comment">%For all questions, bugs and suggestions please mail</span>
0016 <span class="comment">%ivan.oseledets@gmail.com</span>
0017 <span class="comment">%---------------------------</span>
0018 n=a.n;
0019 m=a.m;
0020 att=a.tt;
0021 corea=att.core;
0022 psa=att.ps;
0023 ra=att.r;
0024 d=att.d;
0025 <span class="keyword">if</span> ( nargin &lt;= 5 || isempty(rmax) )
0026    rmax=1000;
0027 <span class="keyword">end</span>
0028 <span class="keyword">if</span> ( nargin &lt;= 4 || isempty(z) )
0029   <span class="comment">%z=tt_random(n,ndims(x),5);</span>
0030   <span class="comment">%z=tt_tensor(z);</span>
0031   z=x;
0032 <span class="keyword">end</span>
0033 <span class="keyword">if</span> ( nargin &lt;= 3 || isempty(nswp) )
0034   nswp = 40;
0035 <span class="keyword">end</span>
0036 y=z;
0037 
0038 <span class="comment">%Warmup is to orthogonalize Y from right-to-left and compute psi-matrices</span>
0039 <span class="comment">%for Ax</span>
0040 psi=cell(d+1,1); <span class="comment">%Psi-matrices</span>
0041 psi{d+1}=1; psi{1}=1;
0042 
0043 
0044 <span class="comment">%Parameters section</span>
0045 
0046 <span class="comment">%Warmup: right-to-left QR + computation of psi matrices</span>
0047 
0048   corex=x.core;
0049   psx=x.ps;
0050   rx=x.r;
0051 
0052     psy=y.ps;
0053   ry=y.r;
0054   corey=y.core;
0055    pos1=psy(d+1);
0056 cr1=corey(psy(d):psy(d+1)-1);
0057 <span class="keyword">for</span> i=d:-1:2  
0058    cr2=corey(psy(i-1):psy(i)-1);
0059    cr1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr1,[ry(i),n(i)*ry(i+1)]);
0060    cr2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr2,[ry(i-1)*n(i-1),ry(i)]);
0061    cr1=cr1.';
0062    [q,rm]=<a href="../../tt2/@tt_tensor/qr.html" class="code" title="function [tt,rm]=qr(tt,op)">qr</a>(cr1,0); rn=<a href="size.html" class="code" title="function [sz] = size(tt)">size</a>(q,2); rm=rm.';
0063    q=q.'; 
0064    ry(i)=rn;
0065    corey(pos1-ry(i+1)*n(i)*ry(i):pos1-1)=q(:);
0066    <span class="comment">%Convolution is now performed for psi(i) using psi(i+1) and corea, and</span>
0067    <span class="comment">%(new) core q</span>
0068    cra=corea(psa(i):psa(i+1)-1); cra=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cra,[ra(i),n(i),m(i),ra(i+1)]);
0069    cry=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(<a href="conj.html" class="code" title="function [b]=conj(a)">conj</a>(q),[ry(i),n(i),ry(i+1)]);
0070    crx=corex(psx(i):psx(i+1)-1); crx=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(crx,[rx(i),m(i),rx(i+1)]);
0071    pscur=psi{i+1}; pscur=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(pscur,[ra(i+1),rx(i+1),ry(i+1)]); <span class="comment">%ra,rx,ry</span>
0072    <span class="comment">%First, convolve over rx(i+1)</span>
0073    crx=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(crx,[rx(i)*m(i),rx(i+1)]);
0074    pscur=permute(pscur,[2,1,3]); pscur=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(pscur,[rx(i+1),ra(i+1)*ry(i+1)]);
0075    pscur=crx*pscur; <span class="comment">%pscur is now rx(i)*m(i)*ra(i+1)*ry(i+1)</span>
0076    <span class="comment">%Convolve over m(i),ra(i+1),n(i),ry(i+1)</span>
0077     pscur=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(pscur,[rx(i),m(i)*ra(i+1),ry(i+1)]);
0078     pscur=permute(pscur,[1,3,2]); 
0079     pscur=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(pscur,[rx(i)*ry(i+1),m(i)*ra(i+1)]);
0080     cra=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cra,[ra(i)*n(i),m(i)*ra(i+1)]); cra=cra.';  
0081     pscur=pscur*cra; 
0082     <span class="comment">%pscur is now rx(i)*ry(i+1)*ra(i)*n(i), it is left to convolve over</span>
0083     <span class="comment">%n(i)*ry(i+1)</span>
0084     pscur=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(pscur,[rx(i),ry(i+1),ra(i),n(i)]);
0085     pscur=permute(pscur,[3,1,4,2]);
0086     pscur=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(pscur,[rx(i)*ra(i),n(i)*ry(i+1)]);
0087     cry=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cry,[ry(i),n(i)*ry(i+1)]); cry=cry.';
0088     pscur=pscur*cry;
0089     psi{i}=pscur;
0090    <span class="comment">%End of psi-block</span>
0091    pos1=pos1-ry(i+1)*n(i)*ry(i);
0092    cr1=cr2*rm;
0093    
0094 <span class="keyword">end</span>
0095 corey(pos1-ry(2)*n(1)*ry(1):pos1-1)=cr1(:);
0096 pos1=pos1-ry(2)*n(1)*ry(1);
0097 corey=corey(pos1:numel(corey)); <span class="comment">%Truncate unused elements</span>
0098 
0099   
0100 psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
0101 y.core=corey;
0102 y.r=ry;
0103 y.ps=psy;
0104 swp=1;
0105    converged=false;
0106 <span class="keyword">while</span> (swp &lt;= nswp &amp;&amp; ~converged) 
0107   <span class="comment">%left-to-right dmrg sweep</span>
0108 psy=y.ps;
0109   ry=y.r;
0110   corey=y.core;
0111   
0112   pos1=1;
0113   cry_old=corey;
0114   converged=true;
0115   ermax=0;
0116   
0117   <span class="keyword">for</span> i=1:d-1
0118      <span class="comment">%We care for two cores, with number i &amp; number i+1, and use</span>
0119      <span class="comment">%psi(i) and psi(i+2) as a basis; also we will need to recompute</span>
0120      <span class="comment">%psi(i+1)</span>
0121      ps1=psi{i}; ps2=psi{i+2};
0122      cra1=corea(psa(i):psa(i+1)-1); cra2=corea(psa(i+1):psa(i+2)-1);
0123      crx1=corex(psx(i):psx(i+1)-1); crx2=corex(psx(i+1):psx(i+2)-1);
0124      <span class="comment">%our convolution is</span>
0125      <span class="comment">%ps1(ra(i),rx(i),ry(i))*cra1(ra(i),n(i),m(i),ra(i+1))*</span>
0126      <span class="comment">%cra2(ra(i+1),n(i+1),m(i+1),ra(i+2))*</span>
0127      <span class="comment">%*ps2(ra(i+2),rx(i+2),ry(i+2))</span>
0128      <span class="comment">%*crx1(rx(i),m(i),rx(i+1))*cr2x(rx(i+1)*m(i+1)*rx(i+2))</span>
0129      <span class="comment">%Scheme ps1*crx1 over rx(i)</span>
0130      
0131      ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ra(i),rx(i),ry(i)]); 
0132      ps1=permute(ps1,[1,3,2]); ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ra(i)*ry(i),rx(i)]);
0133      crx1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(crx1,[rx(i),m(i)*rx(i+1)]);
0134      ps1=ps1*crx1; <span class="comment">%ps1 is now ra(i)*ry(i)*m(i)*rx(i+1)</span>
0135      <span class="comment">%Now convolve with matrix A over ra(i)*m(i)</span>
0136      ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ra(i),ry(i),m(i),rx(i+1)]);
0137      ps1=permute(ps1,[2,4,1,3]);
0138      ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ry(i)*rx(i+1),ra(i)*m(i)]);
0139      cra1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cra1,[ra(i),n(i),m(i),ra(i+1)]);
0140      cra1=permute(cra1,[1,3,2,4]); 
0141      cra1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cra1,[ra(i)*m(i),n(i)*ra(i+1)]);
0142      ps1=ps1*cra1; <span class="comment">%ps1 is now ry(i)*rx(i+1)*n(i)*ra(i+1)</span>
0143      <span class="comment">%Then the ``same'' convolution is carried over for second pair of</span>
0144      <span class="comment">%cores</span>
0145      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[ra(i+2),rx(i+2),ry(i+2)]);
0146      ps2=permute(ps2,[2,1,3]);
0147      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[rx(i+2),ra(i+2)*ry(i+2)]);
0148      crx2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(crx2,[rx(i+1)*m(i+1),rx(i+2)]);
0149      ps2=crx2*ps2; <span class="comment">%ps2 is now rx*(i+1)*m(i+1)*ra(i+2)*ry(i+2)</span>
0150      <span class="comment">%Convolve over m(i+1)*ra(i+2)</span>
0151      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[rx(i+1),m(i+1)*ra(i+2),ry(i+2)]);
0152      ps2=permute(ps2,[2,1,3]);
0153      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[m(i+1)*ra(i+2),rx(i+1)*ry(i+2)]);
0154      cra2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cra2,[ra(i+1)*n(i+1),m(i+1)*ra(i+2)]);
0155      ps2=cra2*ps2; 
0156      <span class="comment">%ps2 is now ra(i+1)*n(i+1)*rx(i+1)*ry(i+2)</span>
0157      <span class="comment">%Now form superblock by contraction ps1 &amp; ps2</span>
0158      <span class="comment">%over ra(i+1)*rx(i+1)</span>
0159      ps0=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ry(i),rx(i+1),n(i),ra(i+1)]);
0160      ps0=permute(ps0,[1,3,2,4]);
0161      ps0=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps0,[ry(i)*n(i),rx(i+1)*ra(i+1)]);
0162      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[ra(i+1),n(i+1),rx(i+1),ry(i+2)]);
0163      ps2=permute(ps2,[3,1,2,4]); 
0164      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[rx(i+1)*ra(i+1),n(i+1)*ry(i+2)]);
0165      super_core=ps0*ps2; <span class="comment">%super_core is ry(i)*n(i)*n(i+1)*ry(i+2)</span>
0166      <span class="comment">%Compute previous supercore</span>
0167      cr1=corey(pos1:pos1+ry(i)*n(i)*ry(i+1)-1);
0168      cr2=cry_old(psy(i+1):psy(i+2)-1); 
0169      cr1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr1,[ry(i)*n(i),ry(i+1)]);
0170      cr2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cr2,[ry(i+1),n(i+1)*ry(i+2)]);
0171      super_core_old=cr1*cr2; 
0172      er=<a href="norm.html" class="code" title="function [nrm] = norm(t,varargin)">norm</a>(super_core_old(:)-super_core(:))/<a href="norm.html" class="code" title="function [nrm] = norm(t,varargin)">norm</a>(super_core(:));
0173      <span class="keyword">if</span> ( er &gt; eps) 
0174         converged=false;
0175      <span class="keyword">end</span>
0176      ermax=max(er,ermax);
0177      <span class="comment">%fprintf('i=%d er=%3.2e \n',i,er);</span>
0178      [u,s,v]=svd(super_core,<span class="string">'econ'</span>);
0179      s=<a href="diag.html" class="code" title="function [tt]=diag(tm)">diag</a>(s); 
0180      r=<a href="../../tt2/core/my_chop2.html" class="code" title="function [r] = my_chop2(sv,eps)">my_chop2</a>(s,eps/sqrt(d-1)*<a href="norm.html" class="code" title="function [nrm] = norm(t,varargin)">norm</a>(s)); r=min(r,rmax);
0181      u=u(:,1:r); s=s(1:r); v=v(:,1:r); v=v*<a href="diag.html" class="code" title="function [tt]=diag(tm)">diag</a>(s);
0182      ry(i+1)=r;
0183      <span class="comment">%u is ry(i)*n(i)*ry(i+1)</span>
0184      <span class="comment">%core_new(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=u(:);</span>
0185      corey(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=u(:); 
0186      
0187      u=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(u,[ry(i),n(i),ry(i+1)]);
0188 
0189      <span class="comment">%Compute new psi</span>
0190      <span class="comment">%ps1 is ry(i)*rx(i+1)*n(i)*ra(i+1) with u over ry(i)*n(i)</span>
0191      u=<a href="conj.html" class="code" title="function [b]=conj(a)">conj</a>(u);
0192      ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ry(i),rx(i+1),n(i),ra(i+1)]);
0193      ps1=permute(ps1,[4,2,1,3]);
0194      ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ra(i+1)*rx(i+1),ry(i)*n(i)]);
0195      u=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(u,[ry(i)*n(i),ry(i+1)]);
0196      ps1=ps1*u;
0197      psi{i+1}=ps1;
0198      <span class="comment">%Compute (?) new v</span>
0199      pos1=pos1+ry(i)*n(i)*ry(i+1);
0200      v=v';
0201      <span class="comment">%v=reshape(v,[ry(i+1),n(i+1),ry(i+2)]);</span>
0202      corey(pos1:pos1+ry(i+1)*n(i+1)*ry(i+2)-1)=v(:);
0203   <span class="keyword">end</span>
0204   psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
0205 y.core=corey;
0206 y.r=ry;
0207 y.ps=psy;
0208 cry_old=corey;
0209 <span class="comment">%m=numel(cry);</span>
0210 <span class="comment">%cry=[zeros(size(cry)),cry];</span>
0211 <span class="comment">%pos1=pos1+m;</span>
0212 <span class="comment">%cry2=cry_old(psy(i+1):psy(i+2)-1);</span>
0213 corey=corey(psy(d):psy(d+1)-1); <span class="comment">%Start--only two cores are left!</span>
0214 
0215  
0216   <span class="keyword">for</span> i=d-1:-1:1
0217      <span class="comment">%We care for two cores, with number i &amp; number i+1, and use</span>
0218      <span class="comment">%psi(i) and psi(i+2) as a basis; also we will need to recompute</span>
0219      <span class="comment">%psi(i)</span>
0220      ps1=psi{i}; ps2=psi{i+2};
0221      cra1=corea(psa(i):psa(i+1)-1); cra2=corea(psa(i+1):psa(i+2)-1);
0222      crx1=corex(psx(i):psx(i+1)-1); crx2=corex(psx(i+1):psx(i+2)-1);
0223      <span class="comment">%our convolution is</span>
0224      <span class="comment">%ps1(ra(i),rx(i),ry(i))*cra1(ra(i),n(i),m(i),ra(i+1))*</span>
0225      <span class="comment">%cra2(ra(i+1),n(i+1),m(i+1),ra(i+2))*</span>
0226      <span class="comment">%*ps2(ra(i+2),rx(i+2),ry(i+2))</span>
0227      <span class="comment">%*crx1(rx(i),m(i),rx(i+1))*cr2x(rx(i+1)*m(i+1)*rx(i+2))</span>
0228      <span class="comment">%Scheme ps1*crx1 over rx(i)</span>
0229      
0230      ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ra(i),rx(i),ry(i)]); 
0231      ps1=permute(ps1,[1,3,2]); ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ra(i)*ry(i),rx(i)]);
0232      crx1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(crx1,[rx(i),m(i)*rx(i+1)]);
0233      ps1=ps1*crx1; <span class="comment">%ps1 is now ra(i)*ry(i)*m(i)*rx(i+1)</span>
0234      <span class="comment">%Now convolve with matrix A over ra(i)*m(i)</span>
0235      ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ra(i),ry(i),m(i),rx(i+1)]);
0236      ps1=permute(ps1,[2,4,1,3]);
0237      ps1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ry(i)*rx(i+1),ra(i)*m(i)]);
0238      cra1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cra1,[ra(i),n(i),m(i),ra(i+1)]);
0239      cra1=permute(cra1,[1,3,2,4]); 
0240      cra1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cra1,[ra(i)*m(i),n(i)*ra(i+1)]);
0241      ps1=ps1*cra1; <span class="comment">%ps1 is now ry(i)*rx(i+1)*n(i)*ra(i+1)</span>
0242      <span class="comment">%Then the ``same'' convolution is carried over for second pair of</span>
0243      <span class="comment">%cores</span>
0244      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[ra(i+2),rx(i+2),ry(i+2)]);
0245      ps2=permute(ps2,[2,1,3]);
0246      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[rx(i+2),ra(i+2)*ry(i+2)]);
0247      crx2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(crx2,[rx(i+1)*m(i+1),rx(i+2)]);
0248      ps2=crx2*ps2; <span class="comment">%ps2 is now rx*(i+1)*m(i+1)*ra(i+2)*ry(i+2)</span>
0249      <span class="comment">%Convolve over m(i+1)*ra(i+2)</span>
0250      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[rx(i+1),m(i+1)*ra(i+2),ry(i+2)]);
0251      ps2=permute(ps2,[2,1,3]);
0252      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[m(i+1)*ra(i+2),rx(i+1)*ry(i+2)]);
0253      cra2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cra2,[ra(i+1)*n(i+1),m(i+1)*ra(i+2)]);
0254      ps2=cra2*ps2; 
0255      <span class="comment">%ps2 is now ra(i+1)*n(i+1)*rx(i+1)*ry(i+2)</span>
0256      <span class="comment">%Now form superblock by contraction ps1 &amp; ps2</span>
0257      <span class="comment">%over ra(i+1)*rx(i+1)</span>
0258      ps0=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps1,[ry(i),rx(i+1),n(i),ra(i+1)]);
0259      ps0=permute(ps0,[1,3,2,4]);
0260      ps0=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps0,[ry(i)*n(i),rx(i+1)*ra(i+1)]);
0261      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[ra(i+1),n(i+1),rx(i+1),ry(i+2)]);
0262      ps2=permute(ps2,[3,1,2,4]); 
0263      ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[rx(i+1)*ra(i+1),n(i+1)*ry(i+2)]);
0264      super_core=ps0*ps2; <span class="comment">%super_core is ry(i)*n(i)*n(i+1)*ry(i+2)</span>
0265      <span class="comment">%Compute previous supercore</span>
0266      cry1=cry_old(psy(i):psy(i+1)-1);
0267      <span class="comment">%cry2=cry(ry(:ry(i+1)*n(i+1)*ry(i+2));</span>
0268      <span class="comment">%cry1=cry(1:ry(i)*n(i)*ry(i+1));</span>
0269      <span class="comment">%cry2=cry(ry(i)*n(i)*ry(i+1)+1:ry(i)*n(i)*ry(i+1)+ry(i+1)*n(i+1)*ry(i+2));</span>
0270      cry2=corey(1:ry(i+1)*n(i+1)*ry(i+2));
0271      corey(1:ry(i+1)*n(i+1)*ry(i+2))=[];
0272      cry1=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cry1,[ry(i)*n(i),ry(i+1)]);
0273      cry2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(cry2,[ry(i+1),n(i+1)*ry(i+2)]);
0274      super_core_old=cry1*cry2;
0275      er=<a href="norm.html" class="code" title="function [nrm] = norm(t,varargin)">norm</a>(super_core_old(:)-super_core(:))/<a href="norm.html" class="code" title="function [nrm] = norm(t,varargin)">norm</a>(super_core(:));
0276      <span class="keyword">if</span> ( er &gt; eps) 
0277         converged=false;
0278      <span class="keyword">end</span>
0279      er=max(er,ermax);
0280      [u,s,v]=svd(super_core,<span class="string">'econ'</span>);
0281      s=<a href="diag.html" class="code" title="function [tt]=diag(tm)">diag</a>(s); 
0282      r=<a href="../../tt2/core/my_chop2.html" class="code" title="function [r] = my_chop2(sv,eps)">my_chop2</a>(s,eps/sqrt(d-1)*<a href="norm.html" class="code" title="function [nrm] = norm(t,varargin)">norm</a>(s)); r=min(r,rmax);
0283      u=u(:,1:r); s=s(1:r); u=u*<a href="diag.html" class="code" title="function [tt]=diag(tm)">diag</a>(s); v=v(:,1:r);
0284      ry(i+1)=r;
0285      v=v.';
0286        u1=u(:); u1=u1.'; v1=v(:); v1=v1.';
0287      corey=[u1,v1,corey];
0288       <span class="comment">%Compute new psi</span>
0289       <span class="comment">%ps2 is rx(i+1)*ra(i+1)*n(i+1)*ry(i+2)</span>
0290       <span class="comment">%v is ry(i+1)*n(i+1)*ry(i+2)</span>
0291       <span class="comment">%convolve over n(i+1)*ry(i+2);</span>
0292       v = v.';
0293       ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,rx(i+1)*ra(i+1),n(i+1)*ry(i+2));
0294       ps2=ps2*v; <span class="comment">%ps2 is now rx(i+1)*ra(i+1)*ry(i+1)</span>
0295       ps2=<a href="../../tt2/@tt_tensor/reshape.html" class="code" title="function [tt2]=reshape(tt1,sz,eps, rl, rr)">reshape</a>(ps2,[rx(i+1),ra(i+1),ry(i+1)]);
0296       ps2=permute(ps2,[2,1,3]);
0297       psi{i+1}=ps2;
0298   <span class="keyword">end</span>  
0299   swp=swp+1;
0300   psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
0301   <span class="comment">%fprintf('swp=%d er=%3.2e \n',swp,ermax);</span>
0302   y.core=corey;
0303   y.r=ry;
0304   y.ps=psy;
0305 <span class="keyword">end</span>
0306 <span class="comment">%Compute alpha*y0 = Ax - y</span>
0307 <span class="comment">%p=p/tt_dot(y0,y0);</span>
0308 <span class="comment">%keyboard;%</span>
0309 <span class="comment">%end</span>
0310 <span class="comment">%p1=tt_mvdot(core(a),core(x),y0);</span>
0311 <span class="comment">%p2=dot(y,tt_tensor(y0));</span>
0312 <span class="comment">%abs(p1-p2)/abs(p1)</span>
0313 <span class="comment">%keyboard;</span>
0314 
0315 
0316</pre></div>
<hr><address>Generated on Wed 08-Feb-2012 18:20:24 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>